using DataFrames, Plots, Statistics, JSON, Clp
plotlyjs()
include(joinpath(dirname(pwd()),"src/TuLiPa.jl")); # using Dates, JuMP, HiGHS, CSV, Clustering
The WebIO Jupyter extension was not detected. See the WebIO Jupyter integration documentation for more information.
┌ Warning: backend `PlotlyJS` is not installed. └ @ Plots C:\Users\jgrc\.julia\packages\Plots\bMtsB\src\backends.jl:43
We show how you can set up a simple deterministic power market problem with TuLiPa. The problem will minimize the cost of covering the demand with production from different technologies.
The problem contains:
These elements need their own model objects. Each model object (e.g. hydropower production or the water balance) has an id, horizon (with time resolution) and special properties. The properties of the object contain information about which elements must be built and updated in the LP problem, as well as a dataset. The dataset can contain values for different weather scenarios, and different model years (e.g. installed capacity in 2021 or 2025). When we update the problem with a weather scenario (scenario time) and a model year (data time), the model objects will build the LP problem for this input data. Then the problem can be solved.
In this demo we use dummy data, and profiles for wind, solar, demand and inflow from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/
The model objects are assembled from a list of dataelements. The dataelements have a concept, concepttype, unique instancename and a data field. The conceptname and concepttype decide how the dataelement should be converted to a model object, while the data field holds all the data that should be put into the model object.
The data elements are stored in a flat structure where each data element holds its information and possibly references to other dataelements they are linked to (for example can an upper capacity hold references to its parameter (stored in another data element), and references to the variable this upper capacity should apply to).
# Read DataElement from a tuple
function getelement(concept, concrete, instance, pairs...)
d = Dict()
for (k, v) in pairs
d[k] = v
end
DataElement(concept, concrete, instance, d)
end
# Power markets or water balances are represented with a Balance equation
# - They have a commodity which will decide the horizon (time-resolution) of the Balance
function addbalance!(elements, name, commodity)
push!(elements, getelement(BALANCE_CONCEPT, "BaseBalance", name,
(COMMODITY_CONCEPT, commodity)))
# Power Balances needs a slack variable if inelastic wind, solar, or run-of-river is higher than the inelastic demand
if commodity == "Power"
slackname = "SlackVar" * name
# Flows are variables that contribute into Balance equations
push!(elements, getelement(FLOW_CONCEPT, "BaseFlow", slackname))
slackarrowname = "SlackArrow" * name
# Arrows connect Flows and Balances
push!(elements, getelement(ARROW_CONCEPT, "BaseArrow", slackarrowname,
(BALANCE_CONCEPT, name),
(FLOW_CONCEPT, slackname),
(CONVERSION_CONCEPT, 1.0), # Factor to convert Flow into the same commodity as Balance (here Power to Power, so 1)
(DIRECTIONKEY, DIRECTIONOUT))) # Positive or negative contribution
end
end
# Rhsterms contribute to the right hand side of a Balance equation
function addrhsterm!(elements, name, balance, direction)
d = getelement(RHSTERM_CONCEPT, "BaseRHSTerm", name,
(BALANCE_CONCEPT, balance),
(PARAM_CONCEPT, name), # constant or time-series data
(DIRECTIONKEY, direction)) # positive or negative contriution to the balance
push!(elements, d)
end
# DataElements for a Thermal power plant
function addrhsthermal!(elements, name, balance;
cost::Union{Real, Nothing}=nothing, # €/GWh
cap::Union{Real, Nothing}=nothing) # MW
genname = "Gen" * name
costname = "MC" * name
capname = "Cap" * name
powerarrowname = "PowerSupply" * name
# Flows are variables that contribute into Balance equations
push!(elements, getelement(FLOW_CONCEPT, "BaseFlow", genname))
# Arrows connect Flows and Balances
push!(elements, getelement(ARROW_CONCEPT, "BaseArrow", powerarrowname,
(BALANCE_CONCEPT, balance),
(FLOW_CONCEPT, genname),
(CONVERSION_CONCEPT, 1.0), # Factor to convert Flow into the same commodity as Balance (here Power to Power, so 1)
(DIRECTIONKEY, DIRECTIONIN))) # Positive or negative contribution
# Cost is the contribution of the variable into the objective function
costdata = getelement(COST_CONCEPT, "CostTerm", costname,
(WHICHCONCEPT, FLOW_CONCEPT),
(WHICHINSTANCE, genname),
(DIRECTIONKEY, DIRECTIONIN)) # positive or negative contribution
if cost === nothing # Point to a parameter DataElement with time-series data
costdata.value[PARAM_CONCEPT] = costname # assume param with same name
else # If the cost is constant for all scenarios
cost = float(cost)
@assert cost >= 0
costdata.value[PARAM_CONCEPT] = cost
end
push!(elements, costdata)
# Non-negative capacity
capdata = getelement(CAPACITY_CONCEPT, "PositiveCapacity", capname,
(WHICHCONCEPT, FLOW_CONCEPT),
(WHICHINSTANCE, genname),
(PARAM_CONCEPT, capname), # Point to a parameter DataElement with time-series data
(BOUNDKEY, BOUNDUPPER)) # Upper or lower capacity
push!(elements, capdata)
if cap != nothing # If the capacity is constant for all scenarios
cap = float(cap)
@assert cap >= 0
# Parameter that converts the capacity in MW to GWh based on the duration of the horizon periods
push!(elements, getelement(PARAM_CONCEPT,"MWToGWhSeriesParam",capname,
("Level",cap),
("Profile",1.0)))
end
return
end
# DataElements for a simple reservoir hydropower plant
function addhydro!(elements, name, powerbalance;
eneq::Union{Real, Nothing}=nothing, # GWh/Mm3 (or kWh/m3)
releasecap::Union{Real, Nothing}=nothing, # m3/s
storagecap::Union{Real, Nothing}=nothing) # Mm3
hydrobalance = name
releasename = "Release" * name
spillname = "Spill" * name
storagename = "Storage" * name
storagecapname = "Cap" * storagename
releasecapname = "Cap" * releasename
inflowname = "Inflow" * hydrobalance
powerarrowname = "PowerSupply" * releasename
hydroarrowname = "WaterDemand" * releasename
spillarrowname = "Spill" * releasename
# Water balance
addbalance!(elements, hydrobalance, "Hydro")
# Inflow is a contribution to the right-hand-side of the water balance
addrhsterm!(elements, inflowname, hydrobalance, DIRECTIONIN)
# Variables for release and spill
push!(elements, getelement(FLOW_CONCEPT, "BaseFlow", releasename))
push!(elements, getelement(FLOW_CONCEPT, "BaseFlow", spillname))
# Release and spill take water out from the hydro balance.
push!(elements, getelement(ARROW_CONCEPT, "BaseArrow", hydroarrowname,
(BALANCE_CONCEPT, hydrobalance),
(FLOW_CONCEPT, releasename),
(CONVERSION_CONCEPT, 1.0),
(DIRECTIONKEY, DIRECTIONOUT)))
push!(elements, getelement(ARROW_CONCEPT, "BaseArrow", spillarrowname,
(BALANCE_CONCEPT, hydrobalance),
(FLOW_CONCEPT, spillname),
(CONVERSION_CONCEPT, 1.0),
(DIRECTIONKEY, DIRECTIONOUT)))
# Release also contributes to the power market
powerarrowdata = getelement(ARROW_CONCEPT, "BaseArrow", powerarrowname,
(BALANCE_CONCEPT, powerbalance),
(FLOW_CONCEPT, releasename),
(DIRECTIONKEY, DIRECTIONIN))
if eneq === nothing
powerarrowdata.value[CONVERSION_CONCEPT] = powerarrowname
else
eneq = float(eneq)
@assert eneq >= 0
powerarrowdata.value[CONVERSION_CONCEPT] = eneq
end
push!(elements, powerarrowdata)
# Release capacity
releasecapdata = getelement(CAPACITY_CONCEPT, "PositiveCapacity", releasecapname,
(WHICHCONCEPT, FLOW_CONCEPT),
(WHICHINSTANCE, releasename),
(PARAM_CONCEPT, releasecapname), # Point to a parameter DataElement with time-series data
(BOUNDKEY, BOUNDUPPER)) # Upper or lower capacity
push!(elements, releasecapdata)
if releasecap != nothing # If the capacity is constant for all scenarios
releasecap = float(releasecap)
@assert releasecap >= 0
# Parameter that converts the capacity in m3/s to Mm3 based on the duration of the horizon periods
push!(elements, getelement(PARAM_CONCEPT,"M3SToMM3SeriesParam",releasecapname,
("Level",releasecap),
("Profile",1.0)))
end
# Variable for storage
push!(elements, getelement(STORAGE_CONCEPT, "BaseStorage", storagename,
(BALANCE_CONCEPT, hydrobalance)))
# Storage capacity
storagecapdata = getelement(CAPACITY_CONCEPT, "PositiveCapacity", storagecapname,
(WHICHCONCEPT, STORAGE_CONCEPT),
(WHICHINSTANCE, storagename),
(BOUNDKEY, BOUNDUPPER))
if storagecap === nothing
storagecapdata.value[PARAM_CONCEPT] = storagecapname
else
storagecap = float(storagecap)
@assert storagecap >= 0
storagecapdata.value[PARAM_CONCEPT] = storagecap
end
push!(elements, storagecapdata)
return
end
# DataElements for transmission between areas
function addpowertrans!(elements, frombalance, tobalance;
cap::Union{Real, Nothing}=nothing,
eff::Union{Real, Nothing}=nothing)
flowname = frombalance * "->" * tobalance
capname = "Cap" * flowname
fromarrowname = flowname * "From"
toarrowname = flowname * "To"
# Transmission variable
push!(elements, getelement(FLOW_CONCEPT, "BaseFlow", flowname))
# Variable out from one Balance
fromarrowdata = getelement(ARROW_CONCEPT, "BaseArrow", fromarrowname,
(BALANCE_CONCEPT, frombalance),
(FLOW_CONCEPT, flowname),
(CONVERSION_CONCEPT, 1.0),
(DIRECTIONKEY, DIRECTIONOUT))
push!(elements, fromarrowdata)
# Variable in to another Balance
toarrowdata = getelement(ARROW_CONCEPT, "BaseArrow", toarrowname,
(BALANCE_CONCEPT, tobalance),
(FLOW_CONCEPT, flowname),
(DIRECTIONKEY, DIRECTIONIN))
if eff === nothing
toarrowdata.value[CONVERSION_CONCEPT] = powerarrowname
else
@assert 0 < eff <= 1
eff = float(eff)
toarrowdata.value[CONVERSION_CONCEPT] = eff # this could also be modelled as a Loss
end
push!(elements, toarrowdata)
# Transmission capacity
capdata = getelement(CAPACITY_CONCEPT, "PositiveCapacity", capname,
(WHICHCONCEPT, FLOW_CONCEPT),
(WHICHINSTANCE, flowname),
(PARAM_CONCEPT, capname), # Point to a parameter DataElement with time-series data
(BOUNDKEY, BOUNDUPPER)) # Upper or lower capacity
push!(elements, capdata)
if cap != nothing # If the capacity is constant for all scenarios
cap = float(cap)
@assert cap >= 0
# Parameter that converts the capacity in MW to GWh based on the duration of the horizon periods
push!(elements, getelement(PARAM_CONCEPT,"MWToGWhSeriesParam",capname,
("Level",cap),
("Profile",1.0)))
end
end
printdicts(elements) = JSON.print(elements, 2);
# Combine the different parts of the dataset into one list of DataElements
function gettestdataset()
elements = DataElement[]
structure = getteststructure()
elements = vcat(elements, structure)
params = gettestparams()
elements = vcat(elements, params)
constants = gettestconstants()
elements = vcat(elements, constants)
levels = gettestlevels()
elements = vcat(elements, levels)
profiles = gettestprofiles()
elements = vcat(elements, profiles)
return elements
end
# The structure consist of the main model objects and how they are connected together
# We also add some of the parameters (like capacities, conversions and costs) if they are constant for all scenarios
function getteststructure()
structure = DataElement[]
addbalance!(structure, "PowerBalance_NO2", "Power")
addbalance!(structure, "PowerBalance_GER", "Power")
addpowertrans!(structure, "PowerBalance_NO2", "PowerBalance_GER", cap=1400, eff=0.97)
addpowertrans!(structure, "PowerBalance_GER", "PowerBalance_NO2", cap=1400, eff=0.97)
addrhsterm!(structure, "WindNO2", "PowerBalance_NO2", DIRECTIONIN)
addrhsterm!(structure, "RoRNO2", "PowerBalance_NO2", DIRECTIONIN)
addrhsterm!(structure, "DemandNO2", "PowerBalance_NO2", DIRECTIONOUT)
addrhsterm!(structure, "WindGER", "PowerBalance_GER", DIRECTIONIN)
addrhsterm!(structure, "SolarGER", "PowerBalance_GER", DIRECTIONIN)
addrhsterm!(structure, "DemandGER", "PowerBalance_GER", DIRECTIONOUT)
addrhsthermal!(structure, "BioGER", "PowerBalance_GER", cap=5000, cost=50000)
addrhsthermal!(structure, "NucGER", "PowerBalance_GER", cap=5000, cost=5000)
addrhsthermal!(structure, "CoalGER", "PowerBalance_GER", cap=15000)
addrhsthermal!(structure, "GasGER", "PowerBalance_GER", cap=40000)
addhydro!(structure, "ResNO2", "PowerBalance_NO2", eneq=1.3, storagecap=16000, releasecap=6000)
return structure
end
# These parameters are built up with TimeVectors (constants or time-series data), which they point to
function gettestparams()
params = DataElement[]
# MWToGWhSeriesParam have a profile and a level stored in other DataElements
# And convert the value from MW to GWh based on the duration of horizon periods
for name in ["WindNO2", "DemandNO2", "WindGER", "SolarGER", "DemandGER"]
push!(params, DataElement(PARAM_CONCEPT, "MWToGWhSeriesParam", name,
Dict("Profile" => "Profile" * name, "Level" => "Level" * name)))
end
# M3SToMM3SeriesParam have a profile and a level stored in other DataElements
# And convert the value from m3/s to Mm3 based on the duration of horizon periods
for name in ["InflowResNO2", "RoRNO2"]
push!(params, DataElement(PARAM_CONCEPT, "M3SToMM3SeriesParam", name,
Dict("Profile" => "Profile" * name, "Level" => "Level" * name)))
end
# FossilMCParam is calculated from many inputs. Here, most of them are constant, except for the
# gas price level and profile that is dependant on chosen scenario
for fuel in ["Coal", "Gas"]
name = "MC" * fuel * "GER"
push!(params, DataElement(PARAM_CONCEPT, "FossilMCParam", name,
Dict("FuelProfile" => "Profile" * fuel, "FuelLevel" => "Level" * fuel,
"CO2Factor" => "CO2Factor" * fuel, "CO2Profile" => "CO2Profile", "CO2Level" => "CO2Level",
"Efficiency" => "Efficiency" * fuel, "VOC" => "VOC" * fuel)))
end
params
end
# We define the constant TimeVectors which make up FossilMCParam
function gettestconstants()
constants = DataElement[]
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "LevelCoal", Dict("Value" => 30000.0))) # €/GWh
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "ProfileCoal", Dict("Value" => 1.0)))
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "CO2FactorGas", Dict("Value" => 0.18)))
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "CO2FactorCoal", Dict("Value" => 0.36)))
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "VOCGas", Dict("Value" => 2000.0))) # €/GWh
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "VOCCoal", Dict("Value" => 4000.0))) # €/GWh
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "CO2Level", Dict("Value" => 50000.0))) # €/GWh
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "CO2Profile", Dict("Value" => 1.0)))
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "EfficiencyGas", Dict("Value" => 0.5)))
push!(constants, DataElement(TIMEVECTOR_CONCEPT, "ConstantTimeVector", "EfficiencyCoal", Dict("Value" => 0.4)))
constants
end
# We define TimeVectors that represent the level of different parameters in 2021 and 2025.
function gettestlevels()
levels = DataElement[]
push!(levels, DataElement(TIMEINDEX_CONCEPT, "VectorTimeIndex", "DataLevelsTimeIndex",
Dict("Vector" => [getisoyearstart(2021), getisoyearstart(2025)])))
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelWindNO2",
Dict("Vector" => Float64[1500, 1500]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelRoRNO2",
Dict("Vector" => Float64[500, 500]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelDemandNO2",
Dict("Vector" => Float64[5000, 5500]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelDemandGER",
Dict("Vector" => Float64[50000, 55000]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelWindGER",
Dict("Vector" => Float64[30000, 40000]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelSolarGER",
Dict("Vector" => Float64[60000, 80000]))) # MW
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelGas",
Dict("Vector" => Float64[90000, 50000]))) # €/GWh
push!(levels, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "LevelInflowResNO2",
Dict("Vector" => Float64[950, 1000]))) # m3/s
for name in ["LevelWindNO2", "LevelRoRNO2", "LevelDemandNO2", "LevelDemandGER",
"LevelWindGER", "LevelSolarGER", "LevelGas", "LevelInflowResNO2"]
push!(levels, DataElement(TIMEVECTOR_CONCEPT, "InfiniteTimeVector", name,
Dict(TIMEINDEX_CONCEPT => "DataLevelsTimeIndex", TIMEVALUES_CONCEPT => name)))
end
levels
end
# We define TimeVectors that represent the profile of different parameters
function gettestprofiles()
# Profiles from https://www.nve.no/energi/analyser-og-statistikk/vaerdatasett-for-kraftsystemmodellene/
path = "testprofiles_1981_2010.csv"
dfmt = dateformat"yyyy-mm-dd HH:MM:SS"
df = CSV.read(path, DataFrame)
df.Timestamp = DateTime.(df.Timestamp, dfmt)
@assert issorted(df.Timestamp)
start = first(df.Timestamp)
numperiods = length(df.Timestamp)
colnames = [n for n in names(df) if n != "Timestamp"]
matrix = Matrix{Float64}(df[:, colnames])
elements = DataElement[]
colnames = ["Profile" * name for name in colnames]
# The time-series data is stored efficiently in two DataElements
# - RangeTimeIndex is a time series index described by a StepRange
# - BaseTable stores the data in a matrix, where each column is a different profile,
# and each row represent the value at a TimeIndex
push!(elements, DataElement(TIMEINDEX_CONCEPT, "RangeTimeIndex", "ProfilesTimeIndex",
Dict("Start" => start, "Delta" => Hour(1), "Steps" => numperiods)))
push!(elements, DataElement(TABLE_CONCEPT, "BaseTable", "ProfilesTable",
Dict("Matrix" => matrix, "Names" => colnames)))
for name in colnames
# ColumnTimeValues points to a row in BaseTable in another DataElement
push!(elements, DataElement(TIMEVALUES_CONCEPT, "ColumnTimeValues", name,
Dict(TABLE_CONCEPT => "ProfilesTable", "Name" => name)))
# RotatingTimeVector points to a index and values pair in other DataElements
push!(elements, DataElement(TIMEVECTOR_CONCEPT, "RotatingTimeVector", name,
Dict(TIMEVALUES_CONCEPT => name, TIMEINDEX_CONCEPT => "ProfilesTimeIndex")))
end
# We make a simple gas price profile for the variation throughout the year
gasprice = Float64[1.079, 1.037, 0.970, 0.931, 0.915, 0.971, 1.028, 1.070]
datetimes = DateTime[getisoyearstart(1981) + Hour((i-1) * 1113) for i in 1:8]
gasprice = gasprice/mean(gasprice)
push!(elements, DataElement(TIMEINDEX_CONCEPT, "VectorTimeIndex", "IndexProfileGas",
Dict("Vector" => datetimes)))
push!(elements, DataElement(TIMEVALUES_CONCEPT, "VectorTimeValues", "ValuesProfileGas",
Dict("Vector" => gasprice)))
push!(elements, DataElement(TIMEVECTOR_CONCEPT, "OneYearTimeVector", "ProfileGas",
Dict(TIMEVALUES_CONCEPT => "ValuesProfileGas", TIMEINDEX_CONCEPT => "IndexProfileGas")))
return elements
end;
Define some model settings that completes the dataset.
See Demo 2 for the use of a more complex horizon which is built based on the residual load.
function gettestrunsettings()
runsettings = DataElement[]
# Select the weather scenarios 1981-1982 (some of the profiles are 1981-2010, we only want to use two years)
# Used by RotatingTimeVector to reduce the amount of data
scenarioyearstart = 1981
scenarioyearstop = 1983
push!(runsettings, getelement(TIMEPERIOD_CONCEPT, "ScenarioTimePeriod", "ScenarioTimePeriod",
("Start", getisoyearstart(scenarioyearstart)), ("Stop", getisoyearstart(scenarioyearstop))))
# We choose the horizon (time-resolution) of the commodities. We set the duration of the horizons to 3 years
# We want the variables connected to power (daily) to be more detailed than the hydro variables (weekly)
global power_horizon = SequentialHorizon(364*3, Day(1))
global hydro_horizon = SequentialHorizon(52*3, Week(1))
push!(runsettings, getelement(COMMODITY_CONCEPT, "BaseCommodity", "Power",
(HORIZON_CONCEPT, power_horizon)))
push!(runsettings, getelement(COMMODITY_CONCEPT, "BaseCommodity", "Hydro",
(HORIZON_CONCEPT, hydro_horizon)))
# Storages have state-dependant variables that need a boundary condition
# We set the starting storage to be equal to the ending storage, x[0] = x[T] (for horizon where t in 1:T)
push!(runsettings, getelement(BOUNDARYCONDITION_CONCEPT, "StartEqualStop", "StartEqualStop_StorageResNO2",
(WHICHINSTANCE, "StorageResNO2"),
(WHICHCONCEPT, STORAGE_CONCEPT)))
return runsettings
end
gettestrunsettings (generic function with 1 method)
testdata = gettestdataset()
runsettings = gettestrunsettings()
elements = vcat(testdata, runsettings)
@time modelobjects = getmodelobjects(elements);
2.871373 seconds (3.58 M allocations: 193.185 MiB, 99.54% compilation time)
display(modelobjects)
Dict{Any, Any} with 15 entries:
Id("Flow", "SlackVarPowe… => BaseFlow(Id("Flow", "SlackVarPowerBalance_NO2"),…
Id("Flow", "GenCoalGER") => BaseFlow(Id("Flow", "GenCoalGER"), SequentialHor…
Id("Balance", "PowerBala… => BaseBalance(Id("Balance", "PowerBalance_GER"), B…
Id("Flow", "ReleaseResNO… => BaseFlow(Id("Flow", "ReleaseResNO2"), Sequential…
Id("Flow", "PowerBalance… => BaseFlow(Id("Flow", "PowerBalance_GER->PowerBala…
Id("Balance", "ResNO2") => BaseBalance(Id("Balance", "ResNO2"), BaseCommodi…
Id("Flow", "GenBioGER") => BaseFlow(Id("Flow", "GenBioGER"), SequentialHori…
Id("Flow", "GenNucGER") => BaseFlow(Id("Flow", "GenNucGER"), SequentialHori…
Id("Balance", "PowerBala… => BaseBalance(Id("Balance", "PowerBalance_NO2"), B…
Id("Flow", "SpillResNO2") => BaseFlow(Id("Flow", "SpillResNO2"), SequentialHo…
Id("Flow", "SlackVarPowe… => BaseFlow(Id("Flow", "SlackVarPowerBalance_GER"),…
Id("Flow", "GenGasGER") => BaseFlow(Id("Flow", "GenGasGER"), SequentialHori…
Id("BoundaryCondition", … => StartEqualStop(Id("BoundaryCondition", "StartEqu…
Id("Storage", "StorageRe… => BaseStorage(Id("Storage", "StorageResNO2"), Base…
Id("Flow", "PowerBalance… => BaseFlow(Id("Flow", "PowerBalance_NO2->PowerBala…
display("We display the fields in the struct")
display(fieldnames(typeof(modelobjects[Id("Flow", "ReleaseResNO2")])))
display("The Id shows that the struct is of type Flow")
display(getid(modelobjects[Id("Flow", "ReleaseResNO2")]))
display("The Horizon consist of 1092 periods with timedelta of 1 day (total 3 years)")
display(gethorizon(modelobjects[Id("Flow", "ReleaseResNO2")]))
display("The lower bound is LowerZeroCapacity so that the variable is non-negative")
display(getlb(modelobjects[Id("Flow", "ReleaseResNO2")]))
display("The upper bound is 6000 m3/s, which is stored in a M3SToMM3SeriesParam. M3SToMM3SeriesParam converts the upper bound into Mm3 depending on the horizon of the problem")
display(getub(modelobjects[Id("Flow", "ReleaseResNO2")]))
display("There are two arrows. The first says that the release takes water out from the water balance and the second says that it contributes into the power balance")
for arrow in getarrows(modelobjects[Id("Flow", "ReleaseResNO2")])
display(getid(arrow))
# display(arrow)
end
"We display the fields in the struct"
(:id, :horizon, :ub, :lb, :costs, :sumcost, :arrows, :metadata)
"The Id shows that the struct is of type Flow"
Id("Flow", "ReleaseResNO2")
"The Horizon consist of 1092 periods with timedelta of 1 day (total 3 years)"
SequentialHorizon(SequentialPeriods([(1092, Millisecond(86400000))]), nothing)
"The lower bound is LowerZeroCapacity so that the variable is non-negative"
LowerZeroCapacity()
"The upper bound is 6000 m3/s, which is stored in a M3SToMM3SeriesParam. M3SToMM3SeriesParam converts the upper bound into Mm3 depending on the horizon of the problem"
PositiveCapacity(Id("Capacity", "CapReleaseResNO2"), M3SToMM3SeriesParam{ConstantTimeVector, ConstantTimeVector}(ConstantTimeVector(6000.0), ConstantTimeVector(1.0)), true)
"There are two arrows. The first says that the release takes water out from the water balance and the second says that it contributes into the power balance"
Id("Arrow", "WaterDemandReleaseResNO2")
Id("Arrow", "PowerSupplyReleaseResNO2")
# It is stored in one of the rhsterms of PowerBalance_NO2
rhsterm = getrhsterms(modelobjects[Id("Balance", "PowerBalance_GER")])[1]
# We start with the installed capacity in 2023 and the weather scenario 1981
start = TwoTime(getisoyearstart(2023), getisoyearstart(1981))
# We look at the weekly production in GWh for three years from the starting point
horizon = SequentialHorizon(52*3, Week(1))
# Vectors to store results
i = Vector{typeof(getdatatime(start))}(undef, getnumperiods(horizon))
v = zeros(getnumperiods(horizon))
# Query results for each period in the horizon
for t in 1:getnumperiods(horizon)
querystart = getstarttime(horizon, t, start)
querydelta = gettimedelta(horizon, t)
i[t] = getdatatime(querystart)
v[t] = getparamvalue(rhsterm, querystart, querydelta)
end
# Plot
display(plot(i, v, legend=false, size=(900,300), ylabel = "GWh", title = "Wind power production NO2"))
We build a problem object for a minimization problem. The problem object consist of:
The problem is also initialized. Variables and balances are built, and variables are put into balances and the objective function if coefficients are constant regardless of scenario or horizon period.
mymodel = JuMP.Model(Clp.Optimizer)
set_silent(mymodel)
prob = JuMP_Prob(modelobjects, mymodel)
# prob = HiGHS_Prob(modelobjects) # we could also have used HiGHS_Prob, but printing the mathematical model is not supported there
prob.model
A JuMP Model
Minimization problem with:
Variables: 10141
Objective function type: AffExpr
`AffExpr`-in-`MathOptInterface.EqualTo{Float64}`: 2341 constraints
`VariableRef`-in-`MathOptInterface.GreaterThan{Float64}`: 10140 constraints
`VariableRef`-in-`MathOptInterface.LessThan{Float64}`: 7800 constraints
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: Clp
Names registered in the model: BalancePowerBalance_GER, BalancePowerBalance_NO2, BalanceResNO2, BoundaryConditionEqStartEqualStop_StorageResNO2, FlowGenBioGER, FlowGenCoalGER, FlowGenGasGER, FlowGenNucGER, FlowPowerBalance_GER->PowerBalance_NO2, FlowPowerBalance_NO2->PowerBalance_GER, FlowReleaseResNO2, FlowSlackVarPowerBalance_GER, FlowSlackVarPowerBalance_NO2, FlowSpillResNO2, StorageStartStorageResNO2, StorageStorageResNO2
The hydrorbalance has been initialized. But the inflow is not added since the value depends on the chosen scenario and period in the horizon.
display(prob.model[:BalanceResNO2][2])
Variables are now put into balances and the objective function if coefficients are dependant on the scenario or horizon period. RHSTerms like the inflow are also updated, but it is stored in the problem object rather than updated in the JuMP model straight away. When solve! is called, the problem object will check if any of the RHSTerm has been updated. If that is the case, the JuMP model will be updated with the RHSTerms.
t = TwoTime(getisoyearstart(2021), getisoyearstart(1981))
update!(prob, t)
@time solve!(prob)
getobjectivevalue(prob)
0.062265 seconds (130.85 k allocations: 11.905 MiB)
9.61646303403377e10
Now the hydropower balance also includes the inflow, which is 156 Mm3 in the second week of the problem.
display(prob.model[:BalanceResNO2][2])
t = TwoTime(getisoyearstart(2024), getisoyearstart(1982))
update!(prob, t)
solve!(prob)
display(getobjectivevalue(prob))
display(prob.model[:BalancePowerBalance_GER][2])
8.194441022906743e10
The prices are in EUR/MWh. The price decline from 2025 is a result of the new gas price level in 2025.
plot(-dual.(prob.model[:BalancePowerBalance_GER])/1000)
# plot(value.(prob.model[:FlowGenGasGER])) # Or we could get the daily gas production in GWh
numperiods = getnumperiods(power_horizon)
x = [getdatatime(t) + getstartduration(power_horizon, s) for s in 1:numperiods]
y = [-getcondual(prob, Id("Balance","PowerBalance_GER"), s) for s in 1:numperiods]/1000
plot(x,y)
The function runs several scenarios and plot prices, supplies, demands and hydro storages for each scenario. The function only plots results for chosen price areas.
In Demo 2 a more complex version of this function is used, which support more types of objects. This includes AdaptiveHorizons (horizon built based on residual load), Horizons with offset, SegmentedArrow (for hydropower plants with PQ-curves), BaseAggSupplyCurve (for aggregated power plants) and ExogenBalance (exogen price areas)
# Make list of scenarios
function getscenarios(dt; years)
[TwoTime(getisoyearstart(dt), getisoyearstart(yr)) for yr in years]
end
# Run scenarios and plot results from the whole problem
function runscenarios(scenarios, modelobjects)
runscenarios(scenarios, modelobjects, values(modelobjects))
end
# Run scenarios and plot results from chosen model objects (e.g. only one price area)
function runscenarios(scenarios, modelobjects, resultobjects)
numperiods_powerhorizon = getnumperiods(power_horizon)
numperiods_hydrohorizon = getnumperiods(hydro_horizon)
# We use the datatime for plotting results
dt = getdatatime(scenarios[1])
# Time vector for power_horizon
x1 = [dt + getstartduration(power_horizon, t) for t in 1:numperiods_powerhorizon]
# Time vector for hydro_horizon
x2 = [dt + getstartduration(hydro_horizon, t) for t in 1:numperiods_hydrohorizon]
# Problem can be a HiGHS_Prob or a JuMP_Prob
# prob = HiGHS_Prob(collect(values(modelobjects)))
model = Model(HiGHS.Optimizer)
set_silent(model)
prob = JuMP_Prob(collect(values(modelobjects)), model)
# Order result objects into lists
powerbalances = []
rhsterms = []
plants = []
plantbalances = []
plantarrows = Dict()
demands= []
demandbalances = []
demandarrows = Dict()
hydrostorages = []
for obj in resultobjects
# Powerbalances
if obj isa BaseBalance
if getinstancename(getid(getcommodity(obj))) == "Power"
push!(powerbalances, getid(obj))
for rhsterm in getrhsterms(obj)
push!(rhsterms,getid(rhsterm))
end
end
end
# Hydrostorages
if obj isa BaseStorage
if getinstancename(getid(getcommodity(getbalance(obj)))) == "Hydro"
push!(hydrostorages,getid(obj))
end
end
# Supply and demands
if obj isa BaseFlow
# The type of supply or demand can be found based on the arrows
arrows = getarrows(obj)
# Simple supplies and demands
powerarrowbool = [getid(getcommodity(getbalance(arrow))) == Id("Commodity", "Power") for arrow in arrows]
powerarrows = arrows[powerarrowbool]
if sum(powerarrowbool) == 1
if isingoing(powerarrows[1])
push!(plants,getid(obj))
push!(plantbalances,getid(getbalance(powerarrows[1])))
elseif !isingoing(powerarrows[1])
push!(demands,getid(obj))
push!(demandbalances,getid(getbalance(powerarrows[1])))
end
end
# Transmissions
if sum(powerarrowbool) == 2
for arrow in arrows
balance = getbalance(arrow)
if getid(getcommodity(balance)) == Id("Commodity", "Power")
if isingoing(arrow) && (balance in resultobjects)
push!(plants,getid(obj))
push!(plantbalances,getid(balance))
elseif !isingoing(arrow) && (balance in resultobjects)
push!(demands,getid(obj))
push!(demandbalances,getid(balance))
end
end
end
end
end
end
# Matrices to store results per time period, scenario and object
prices = zeros(numperiods_powerhorizon, length(scenarios), length(powerbalances))
rhstermvalues = zeros(numperiods_powerhorizon, length(scenarios), length(rhsterms))
production = zeros(numperiods_powerhorizon, length(scenarios), length(plants))
consumption = zeros(numperiods_powerhorizon, length(scenarios), length(demands))
hydrolevels = zeros(numperiods_hydrohorizon, length(scenarios), length(hydrostorages))
# Update and solve scenarios, and collect results
for (s, t) in enumerate(scenarios)
update!(prob, t)
solve!(prob)
println("Objective value in scenario $(s): ", getobjectivevalue(prob))
# Collect results for power production/demand/transmission in GW
for j in 1:numperiods_powerhorizon
# Timefactor transform results from GWh to GW/h regardless of horizon period durations
timefactor = getduration(gettimedelta(power_horizon, j))/Millisecond(3600000)
# For powerbalances collect prices and rhsterms (like inelastic demand, wind, solar and RoR)
for i in 1:length(powerbalances)
prices[j, s, i] = -getcondual(prob, powerbalances[i], j)/1000 # from €/GWh to €/MWh
for k in 1:length(rhsterms)
if hasrhsterm(prob, powerbalances[i], rhsterms[k], j)
rhstermvalues[j, s, k] = getrhsterm(prob, powerbalances[i], rhsterms[k], j)/timefactor
end
end
end
# Collect production of all plants
for i in 1:length(plants) # TODO: Balance and variable can have different horizons
production[j, s, i] = getvarvalue(prob, plants[i], j)*abs(getconcoeff(prob, plantbalances[i], plants[i], j, j))/timefactor
end
# Collect demand of all demands
for i in 1:length(demands) # TODO: Balance and variable can have different horizons
consumption[j, s, i] = getvarvalue(prob, demands[i], j)*abs(getconcoeff(prob, demandbalances[i], demands[i], j, j))/timefactor
end
end
# Collect hydro storage levels
for j in 1:numperiods_hydrohorizon
for i in 1:length(hydrostorages)
hydrolevels[j, s, i] = getvarvalue(prob, hydrostorages[i], j)/1000 # Gm3 TODO: convert to TWh with global energy equivalents of each storage
end
end
end
# Only keep rhsterms that have at least one value (TODO: Do the same for sypply and demands)
rhstermtotals = dropdims(sum(rhstermvalues,dims=(1,2)),dims=(1,2))
rhstermsupplyidx = []
rhstermdemandidx = []
for k in 1:length(rhsterms)
if rhstermtotals[k] > 0
push!(rhstermsupplyidx, k)
elseif rhstermtotals[k] < 0
push!(rhstermdemandidx, k)
end
end
# Put rhsterms together with supplies and demands
rhstermsupplyvalues = rhstermvalues[:,:,rhstermsupplyidx]
rhstermdemandvalues = rhstermvalues[:,:,rhstermdemandidx]*-1
rhstermsupplynames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermsupplyidx]]
rhstermdemandnames = [getinstancename(rhsterm) for rhsterm in rhsterms[rhstermdemandidx]]
supplynames = [[getinstancename(plant) for plant in plants];rhstermsupplynames]
supplyvalues = reshape([vcat(production...);vcat(rhstermsupplyvalues...)],(numperiods_powerhorizon,length(scenarios),length(supplynames)))
demandnames = [[getinstancename(demand) for demand in demands];rhstermdemandnames]
demandvalues = reshape([vcat(consumption...);vcat(rhstermdemandvalues...)],(numperiods_powerhorizon,length(scenarios),length(demandnames)))
# Prepare for plotting results
hydronames = [getinstancename(hydro) for hydro in hydrostorages]
powerbalancenames = [split(getinstancename(powerbalance), "PowerBalance_")[2] for powerbalance in powerbalances]
# Plot results for each scenario
for (s, t) in enumerate(scenarios)
scenyear = string(getisoyear(getscenariotime(t)))
datayear = string(getisoyear(getdatatime(t)))
# Plot prices
display(plot(x1, sum(prices[:,s,:],dims=3), labels=reshape(powerbalancenames,1,length(powerbalancenames)), size=(800,500), title="Prices for scenario that starts in " * datayear * " and weather scenario " * scenyear, ylabel="€/MWh"))
# Plot supplies and demands
supplychart = areaplot(x1, sum(supplyvalues[:,s,:],dims=3),labels=reshape(supplynames,1,length(supplynames)),title="Supply", ylabel = "GWh/h")
demandchart = areaplot(x1, sum(demandvalues[:,s,:],dims=3),labels=reshape(demandnames,1,length(demandnames)),title="Demand", ylabel = "GWh/h")
display(plot([supplychart,demandchart]...,layout=(1,2),size=(1600,500)))
# Plot storages (only TWh because of our input data)
display(areaplot(x2, sum(hydrolevels[:,s,:],dims=3),labels=reshape(hydronames,1,length(hydronames)),size=(800,500),title="Reservoir levels", ylabel = "Gm3")) #
# Plot list of yearly mean production and demand for each supply/demand
meandemand = dropdims(mean(demandvalues[:,s,:],dims=1),dims=1)
meanproduction = dropdims(mean(supplyvalues[:,s,:],dims=1),dims=1)
supplydf = sort(DataFrame(Supplyname = supplynames, Yearly_supply_TWh = meanproduction*8.76),[:Yearly_supply_TWh], rev = true)
demanddf = sort(DataFrame(Demandname = demandnames, Yearly_demand_TWh = meandemand*8.76),[:Yearly_demand_TWh], rev = true)
supplydf[!,:ID] = collect(1:length(supplynames))
demanddf[!,:ID] = collect(1:length(demandnames))
joineddf = select!(outerjoin(supplydf,demanddf;on=:ID),Not(:ID))
show(joineddf,allcols=true, allrows=true)
# Check that total supply equals total demand
show(combine(joineddf, [:Yearly_supply_TWh, :Yearly_demand_TWh] .=> sum∘skipmissing))
end
end;
We choose scenarios that start in the model year 2024 and weather scenarios that start in 1981/1982/1983
We get expected results for the deterministic market problem
scenarios = getscenarios(2024; years=1981:1983)
runscenarios(scenarios, modelobjects)
# alternative for only seeing results from one price area
# resultobjects = getpowerobjects(modelobjects,["NO2"])
# runscenarios(scenarios, modelobjects, resultobjects)
Objective value in scenario 1: 8.264928001510184e10 Objective value in scenario 2: 8.19444102290675e10 Objective value in scenario 3: 8.303707657019687e10
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 131.058 DemandGER 468.054 2 │ WindGER 111.944 DemandNO2 48.4256 3 │ SolarGER 70.6705 PowerBalance_NO2->PowerBalance_G… 8.62053 4 │ GenCoalGER 60.4118 SlackVarPowerBalance_NO2 1.80325 5 │ GenNucGER 43.5204 PowerBalance_GER->PowerBalance_N… 0.984684 6 │ GenBioGER 43.1119 SlackVarPowerBalance_GER 0.0400589 7 │ ReleaseResNO2 36.954 missing missing 8 │ RoRNO2 15.9198 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.36191 missing missing 10 │ WindNO2 5.02039 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.955143 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 527.929 527.929
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 129.634 DemandGER 467.293 2 │ WindGER 109.393 DemandNO2 47.626 3 │ SolarGER 73.102 PowerBalance_NO2->PowerBalance_G… 9.10136 4 │ GenCoalGER 60.1686 SlackVarPowerBalance_NO2 1.48872 5 │ GenNucGER 43.5209 PowerBalance_GER->PowerBalance_N… 0.412345 6 │ GenBioGER 43.1053 SlackVarPowerBalance_GER 0.0463198 7 │ ReleaseResNO2 37.2763 missing missing 8 │ RoRNO2 15.4017 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.82832 missing missing 10 │ WindNO2 5.13809 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.399974 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 525.968 525.968
11×4 DataFrame Row │ Supplyname Yearly_supply_TWh Demandname Yearly_demand_TWh │ String? Float64? String? Float64? ─────┼──────────────────────────────────────────────────────────────────────────────────────────────────────────── 1 │ GenGasGER 131.754 DemandGER 468.242 2 │ WindGER 111.142 DemandNO2 48.4851 3 │ SolarGER 70.6877 PowerBalance_NO2->PowerBalance_G… 8.56625 4 │ GenCoalGER 60.7569 SlackVarPowerBalance_NO2 1.79717 5 │ GenNucGER 43.5204 PowerBalance_GER->PowerBalance_N… 0.994443 6 │ GenBioGER 43.1119 SlackVarPowerBalance_GER 0.0463198 7 │ ReleaseResNO2 36.9577 missing missing 8 │ RoRNO2 15.9212 missing missing 9 │ PowerBalance_NO2->PowerBalance_G… 8.30926 missing missing 10 │ WindNO2 5.005 missing missing 11 │ PowerBalance_GER->PowerBalance_N… 0.96461 missing missing 1×2 DataFrame Row │ Yearly_supply_TWh_sum_skipmissing Yearly_demand_TWh_sum_skipmissing │ Float64 Float64 ─────┼────────────────────────────────────────────────────────────────────── 1 │ 528.131 528.131